#pragma once
#include "elisa/elisa.h"
#include <iostream>
#include <algorithm>
#include <string>
#include <vector>
#include <gaol/gaol>

#include "Modeler.h"

#ifdef ELISA_PERFORMER
#include "elisa_chrono.h"
#endif

using namespace std;
using namespace elisa;

/* Problem: Chemical equilibrium of hydrocarbon combustion

@Article{MeintjesMorgan1992,
  author    = {K. Meintjes and A. P. Morgan},
  title     = {Chemical equilibrium systems as numerical test problems},
  year      = {1990},
  journal   = {ACM TOMS},
  volume    = {16},
  number    = {2},
  pages     = {143--151}
}

Although the total degree equals 108, there are only 4 real and
12 complex solutions and an infinite number of solutions at infinity.
*/
class Chemistry2 : public Modeler
{
private:
	Model _model;

	RealConstant Forty;
	RealConstant R;
	RealConstant R5;
	RealConstant R6;
	RealConstant R7;
	RealConstant R8;
	RealConstant R9;
	RealConstant R10;
	RealVar x1;
	RealVar x2;
	RealVar x3;
	RealVar x4;
	RealVar x5;
	RealConstraint c1;
	RealConstraint c2;
	RealConstraint c3;
	RealConstraint c4;
	RealConstraint c5;

public:
	Chemistry2()
		: R(10),
		  R5(0.193),
		  R6(0.0004106), // 0.002597/sqrt(40)
		  R7(0.0005452), // 0.003448/sqrt(40)
		  R8(0.0000004498),
		  R9(0.00003407), // 0.0002155/sqrt(40)
		  R10(0.0000009615),
		  x1(Interval(0,1.0e8), "x1"),
		  x2(Interval(0,1.0e8), "x2"),
		  x3(Interval(0,1.0e8), "x3"),
		  x4(Interval(0,1.0e8), "x4"),
		  x5(Interval(0,1.0e8), "x5"),
		  c1(3*x5 == x1*(x2 + 1)),
		  c2(x3*(x2*(2*x3+R7) + 2*R5*x3 + R6) == 8*x5),
		  c3(x4*(R9*x2 + 2*x4) == 4*R*x5),
		  c4(x2*(2*x1 + x3*(x3+R7) + R8 + 2*R10*x2 + R9*x4) + x1 == R*x5),
		  c5(x2*(x1 + R10*x2 + x3*(x3+R7) + R8 + R9*x4) + x1 + x3*(R5*x3 + R6) + pow(x4,2) == 1),
		  _model()
	{
	}

	elisa::Model & model()
	{
		_model = c1 && c2 && c3 && c4 && c5;
		return _model;
	}

	void print()
	{
		cout << '\t' << x1.getName() << " = " << x1.domain() << endl;
		cout << '\t' << x2.getName() << " = " << x2.domain() << endl;
		cout << '\t' << x3.getName() << " = " << x3.domain() << endl;
		cout << '\t' << x4.getName() << " = " << x4.domain() << endl;
		cout << '\t' << x5.getName() << " = " << x5.domain() << endl;
	}

	RealVar & varX()
	{
		return x1;
	}
	RealVar & varY()
	{
		return x2;
	}
};
